Crop Productivity Analysis#
Data#
The following datasets are utilized in this analysis for calculating and mapping crop productivity over the past years:
Dynamic World Dataset:
Source: Dynamic World - Google and the World Resources Institute (WRI)
Description: The Dynamic World dataset provides a near real-time, high-resolution (10-meter) global land cover classification. It is derived from Sentinel-2 imagery and utilizes machine learning models to classify land cover into nine distinct classes, including water, trees, grass, crops, built areas, bare ground, shrubs, flooded vegetation, and snow/ice. The dataset offers data with minimal latency, enabling near-immediate analysis and decision-making.
Spatial Resolution: 10 meters.
Temporal Coverage: Data is available since mid-2015, updated continuously as Sentinel-2 imagery becomes available capturing near Real-time.
MODIS Dataset:
Source: NASA’s Moderate Resolution Imaging Spectroradiometer MODIS on Terra and Aqua satellites.
Description: The MODIS dataset provides a wide range of data products, including land surface temperature, vegetation indices, and land cover classifications. It is widely used for monitoring and modeling land surface processes.
Spatial Resolution: 250 meters.
Temporal Coverage: Data is available from 2000 to the present, with daily to 16-day composite products.
GDP and Agricultural Data:
Source: Official government statistics and international databases.
Description: This dataset includes regional GDP figures, agricultural output, export data, and crop production statistics. It is used to analyze the economic impact of agricultural productivity.
Spatial Resolution: Varies by dataset; typically available at national and subnational levels.
Temporal Coverage: Varies by dataset; typically available annually or quarterly.
Administrative Boundaries:
Source: MIMU (Myanmar Information Management Unit)
Description: This dataset provides the administrative boundaries for Myanmar at various levels (national, state/region, district, township). It is used for spatial aggregation and analysis of crop productivity and economic data.
Spatial Resolution: Varies by administrative level.
Assumptions and Methodological Notes#
This analysis is based on the following assumptions and methodological considerations:
Temporal Definitions and Aggregations#
Fiscal Year Definition: The first quarter (Q1) of each fiscal year starts in April and ends in March of the following year. This aligns with Myanmar’s fiscal calendar.
Quarterly Aggregations: Monthly data are aggregated into quarters (Q1: April-June, Q2: July-September, Q3: October-December, Q4: January-March).
Methodological Assumptions#
Spatial Aggregation: Regional GDP values are based on single-year reference data and are assumed to maintain consistent spatial distributions across all years in the analysis period.
EVI Processing:
EVI values are median-aggregated within administrative boundaries to reduce noise
Time series preprocessing follows TIMESAT methodology: outlier removal, interpolation, and Savitzky-Golay smoothing
Seasonality parameters (start of season, middle of season, end of season) are extracted using a 20% threshold of maximum EVI
Cropland Identification: Pixels are classified as cropland based on Dynamic World’s machine learning classification. Classification accuracy varies by region and land cover type, with potential confusion between cropland and similar land covers.
def find_project_root(marker_file: str = "pyproject.toml"):
for parent in Path().cwd().parents:
if (parent / marker_file).exists():
return parent
raise FileNotFoundError(f"No project root found with marker file: {marker_file}")
PROJECT_ROOT = Path().cwd().parent.parent
DATA_PATH = PROJECT_ROOT / "data"
EVI_PATH = DATA_PATH / "EVI and Crop Land" / "EVI 2010-2025"
CROPLAND_PATH = DATA_PATH / "EVI and Crop Land" / "Crop Land"
BOUNDARIES_PATH = DATA_PATH / "Shapefiles"
Crop area statistics#
| Name | PCODE | Crop Area (ha) in 2015 | Crop Area (ha) in 2025 | % Change in Crop Area (2015-2025) | Absolute % Change in Crop Area (2015-2025) | |
|---|---|---|---|---|---|---|
| 0 | Sagaing | MMR005 | 897012.27 | 1355659.42 | 51.13 | 51.13 |
| 1 | Ayeyarwady | MMR017 | 477942.32 | 1185667.23 | 148.08 | 148.08 |
| 2 | Mandalay | MMR010 | 802809.95 | 919394.40 | 14.52 | 14.52 |
| 3 | Magway | MMR009 | 739612.49 | 832402.33 | 12.55 | 12.55 |
| 4 | Bago (East) | MMR007 | 3294.37 | 612430.37 | 18490.21 | 18490.21 |
| 5 | Yangon | MMR013 | 1962.76 | 469040.30 | 23796.98 | 23796.98 |
| 6 | Bago (West) | MMR008 | 309351.34 | 442237.44 | 42.96 | 42.96 |
| 7 | Shan (South) | MMR014 | 4505.87 | 300314.13 | 6564.95 | 6564.95 |
| 8 | Mon | MMR011 | 7608.89 | 178099.30 | 2240.67 | 2240.67 |
| 9 | Kachin | MMR001 | 99817.44 | 164590.06 | 64.89 | 64.89 |
| 10 | Shan (North) | MMR015 | 34504.25 | 154754.92 | 348.51 | 348.51 |
| 11 | Rakhine | MMR012 | 178160.82 | 142033.64 | -20.28 | 20.28 |
| 12 | Nay Pyi Taw | MMR018 | 85037.49 | 117129.11 | 37.74 | 37.74 |
| 13 | Kayin | MMR003 | 11609.83 | 82406.40 | 609.80 | 609.80 |
| 14 | Shan (East) | MMR016 | 19354.78 | 36171.34 | 86.89 | 86.89 |
| 15 | Tanintharyi | MMR006 | 11120.88 | 12136.31 | 9.13 | 9.13 |
| 16 | Kayah | MMR002 | 0.00 | 10242.61 | 0.00 | 0.00 |
| 17 | Chin | MMR004 | 6177.81 | 10230.50 | 65.60 | 65.60 |
Crop Seasonality#
Using this time series dataset of EVI images, we apply several pre-processing steps to extract critical phenological parameters: start of season (SOS), middle of season (MOS), end of season (EOS), length of season (LOS), etc. This workflow is heavily inspired by the TIMESAT software.
Pre-processing steps
Remove outliers from dataset on per-pixel basis using median method: outlier if median from a moving window < or > standard deviation of time-series times 2.
Interpolate missing values linearly
Smooth data on per-pixel basis (using Savitsky Golay filter, window length of 3, and polyorder of 1)
Phenology Process
We then extract crop seasonality metrics using the seasonal amplitude method from the phenolopy package.
The chart below shows the result of this process for a single crop pixel. The blue dots represent the raw EVI values, the black line represents the processed EVI values, and the dotted lines represent season parameters extracted for that pixel: start of season, peak of season, and end of season.
Based on the phenology process, we identified the seasonality to start in July and end in February with the peak being in October. This can vary with geographic region and crop type as well, however, that has not been taken into consideration in this version.
The following figure shows the median EVI in Myanmar from 2010 to 2025 during the crop growing season (June to February). The shaded area represents the interquartile range (IQR) of EVI values across all pixels, indicating the variability in vegetation health during this period.
Trends in EVI#
Admin level 1#
The figure below shows the trends of median EVI from 2010-2025 on admin level 1.
The figure below shows a choropleth maps of EVI for each admin level 1 from 2010-2025
The figure below shows the percent change of EVI compared to previous year from 2010-2025 on admin level 1.
Admin level 2#
The figure below shows a choropleth maps of EVI for each admin level 2 from 2010-2025
Admin level 3#
The figure below shows a choropleth maps of EVI for each admin level 3 from 2010-2025